import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.offline as pyo
import plotly.graph_objects as go
import matplotlib.dates as mdates
from plotly.subplots import make_subplots
import plotly.io as pio
pio.renderers.default = 'notebook_connected'
from datetime import timedelta

df = pd.read_csv('simulation_results_finished.csv')
df['Date'] = pd.to_datetime(df['Date'])  # Convert 'Date' column to datetime
df.set_index('Date', inplace=True)
outdoor_temp_series = df['Dry-bulb temperature (°C)']
outdoor_temp_series = outdoor_temp_series.to_frame().reset_index()

# Convert time to seconds for simulation alignment instead of half-hourly intervals
outdoor_temp_series["time_seconds"] = (outdoor_temp_series["Date"] - outdoor_temp_series["Date"].iloc[0]).dt.total_seconds()

# Define constants
heater_power = 2000  # Heater power in watts (W)
room_volume = 180*2.6     # Room volume in cubic meters (m³)
air_density = 1.225  # Air density in kg/m³
specific_heat_air = 1005  # Specific heat of air in J/(kg·K)
heat_loss_rate = 30  # Heat loss rate to environment in W/K (depends on insulation)
initial_indoor_temp = 20  # Initial indoor temperature in °C
thermostat_temp = 20  # Thermostat set temperature in °C
hysteresis = 2.0  # Thermostat hysteresis in °C
simulation_time = 365 * 24 * 3600  # Simulation time for one year (in seconds)
time_step = 1  # Time step in seconds
air_mass_flow_rate = 0.1  # Air mass flow rate in kg/s

# Derived constants
air_mass = room_volume * air_density  # Total air mass in the room (kg)

# Simulation variables
time = np.arange(0, simulation_time, time_step)  # Time array for one year
indoor_temp = np.zeros_like(time, dtype=float)  # indoor temperature array
heater_status = np.zeros_like(time, dtype=int)  # Heater status array (1=ON, 0=OFF)
heater_temp = np.zeros_like(time, dtype=float)  # Heater temperature array
heat_loss = np.zeros_like(time, dtype=float)  # Heat loss array


# Interpolate outdoor temperature for simulation time steps
outdoor_temp = np.interp(time, outdoor_temp_series["time_seconds"], outdoor_temp_series['Dry-bulb temperature (°C)'])

# Initialize indoor temperature
indoor_temp[0] = initial_indoor_temp

# Simulation loop
for t in range(1, len(time)):
    # Calculate heat losses
    heat_loss[t-1] = heat_loss_rate * (indoor_temp[t-1] - outdoor_temp[t-1])
    
    # Determine heater status based on thermostat control
    if indoor_temp[t-1] < thermostat_temp - hysteresis / 2:
        heater_status[t] = 1  # Heater ON
    elif indoor_temp[t-1] > thermostat_temp + hysteresis / 2:
        heater_status[t] = 0  # Heater OFF
    else:
        heater_status[t] = heater_status[t-1]  # Maintain previous state

    # Calculate heat gain from heater
    heat_gain = heater_status[t] * heater_power
    
    # Calculate heater temperature
    heater_temp[t-1] = indoor_temp[t-1] + heat_gain / (air_mass_flow_rate * specific_heat_air)

    # Calculate net heat change
    net_heat = heat_gain - heat_loss[t-1]

    # Update indoor temperature
    temp_change = (net_heat * time_step) / (air_mass * specific_heat_air)
    indoor_temp[t] = indoor_temp[t-1] + temp_change


# Save the results in a DataFrame
sim_results = pd.DataFrame({
    'Time': time,
    'Outdoor Temperature (°C)': outdoor_temp,
    'Indoor Temperature (°C)': indoor_temp,
    'Heater Temperature (°C)': heater_temp,
    'Heat Gain Power (W)': heater_status * heater_power,
    'Heat Gain Energy (kWh)': (heater_status * heater_power * time_step) / 3.6e6,
    'Heat Pump Status': heater_status,
    'Heat Loss Power (W)': heat_loss,
    'Heat Loss Energy (kWh)': (heat_loss * time_step) / 3.6e6,
})

sim_results.set_index('Time', inplace=True)
sim_results.index = pd.to_datetime(sim_results.index, unit='s', origin='unix').strftime('%m-%d %H:%M:%S')
sim_results.to_csv('../Rule based control - heat pumps/dataset/resistive_heater.csv', index=True)

df = pd.read_csv('../Rule based control - heat pumps/dataset/resistive_heater.csv')
df.set_index('Time', inplace=True)
df.head()

# Calculate heater runtime and energy usage
heater_runtime = np.sum(df['Heat Pump Status']) * time_step / 3600  # Hours
energy_usage = np.sum(df['Heat Gain Energy (kWh)']) # kWh

# Display results
print(f"Heater runtime: {heater_runtime:.2f} hours")
print(f"Total energy usage: {energy_usage:.2f} kWh")

Heater runtime: 1562.87 hours
Total energy usage: 3125.74 kWh

# Function to filter data for a specific day
def filter_data_for_day(df, day):
    start_time = pd.to_datetime(day)
    end_time = start_time + pd.Timedelta(days=1)
    return df[(df.index >= start_time.strftime('%m-%d %H:%M:%S')) & (df.index < end_time.strftime('%m-%d %H:%M:%S'))]

# Choose a specific day to display
specific_day = '2023-08-20 00:00:00'  # Change this to the desired date format 'YYYY-MM-DD HH:MM:SS'

daily_profile = filter_data_for_day(df, specific_day)

# Create figure with secondary y-axis
fig = make_subplots(specs=[[{"secondary_y": True}]])

# Add traces
fig.add_trace(
    go.Scatter(x=daily_profile.index, y=daily_profile['Indoor Temperature (°C)'], name='Indoor Temperature (°C)', line=dict(color='blue')),
    secondary_y=False,
)

fig.add_trace(
    go.Scatter(x=daily_profile.index, y=daily_profile['Outdoor Temperature (°C)'], name='Outdoor Temperature (°C)', line=dict(color='orange')),
    secondary_y=False,
)

fig.add_trace(
    go.Scatter(x=daily_profile.index, y=daily_profile['Heat Gain Power (W)'], name='Heat Gain Power (W)', line=dict(color='green', dash='dot')),
    secondary_y=True,
)

# Add thermostat bounds
fig.add_trace(go.Scatter(x=daily_profile.index, y=[thermostat_temp]*len(daily_profile.index), line=dict(color='red', dash='dash'), name="Thermostat Temperature (°C)"), secondary_y=False)
fig.add_trace(go.Scatter(x=daily_profile.index, y=[thermostat_temp - hysteresis / 2]*len(daily_profile.index), line=dict(color='grey', dash='dash'), name="Thermostat Lower Limit (°C)"), secondary_y=False)
fig.add_trace(go.Scatter(x=daily_profile.index, y=[thermostat_temp + hysteresis / 2]*len(daily_profile.index), line=dict(color='grey', dash='dash'), name="Thermostat Upper Limit (°C)"), secondary_y=False)


# Add figure title and axis labels
fig.update_layout(
    title_text=f"Indoor temperature, outdoor temperature, and heat gain power from {pd.to_datetime(specific_day).strftime('%m-%d %H:%M:%S')} to {(pd.to_datetime(specific_day) + timedelta(days=1)).strftime('%m-%d %H:%M:%S')}",
    width=1000,  # Set the width of the figure
    height=700,   # Set the height of the figure
    legend=dict(
        orientation="h",  # Horizontal orientation
        x=0.5,  # x position of the legend
        y=-0.4,  # y position of the legend
        xanchor="center",  # Anchor the legend at the center
        yanchor="bottom",  # Anchor the legend at the bottom
        bgcolor='rgba(255, 255, 255, 0.5)',  # Background color of the legend
        bordercolor='Black',  # Border color of the legend
        borderwidth=1,  # Border width of the legend
        itemwidth=30,  # Width of each legend item
        traceorder="normal"  # Order of legend items
    ),
    xaxis=dict(
        tickformat='%H:%M',  # Format ticks as hours and minutes
        dtick=3600,  # Set tick step to 1 hour (3600000 milliseconds)
        tickangle=45  # Rotate x-axis labels for better readability
    ),
)
# Set x-axis title
fig.update_xaxes(title_text="Time (s)")

# Update y-axes
fig.update_yaxes(title_text='Temperature (°C)', secondary_y=False)
fig.update_yaxes(title_text='Heat Gain Power (W)', secondary_y=True, tickfont=dict(color='green'), titlefont=dict(color='green'))

# Show plot
fig.show()

	Outdoor Temperature (°C)	Indoor Temperature (°C)	Heater Temperature (°C)	Heat Gain Power (W)	Heat Gain Energy (kWh)	Heat Pump Status	Heat Loss Power (W)	Heat Loss Energy (kWh)
Time
01-01 00:00:00	5.5	20.000000	20.000000	0	0.0	0	435.000000	0.000121
01-01 00:00:01	5.5	19.999245	19.999245	0	0.0	0	434.977350	0.000121
01-01 00:00:02	5.5	19.998490	19.998490	0	0.0	0	434.954702	0.000121
01-01 00:00:03	5.5	19.997735	19.997735	0	0.0	0	434.932054	0.000121
01-01 00:00:04	5.5	19.996980	19.996980	0	0.0	0	434.909408	0.000121

Resistive heater hysteresis control¶

References¶

Visualisation¶